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A  ID  model  of  a  plasma  actuator  has  been  developed.  The  plasma  actuator  is  modeled  as 
a  dielectric  barrier  discharge  (DBD)  and  is  based  on  the  numerical  solution  of  the  electron 
and  ion  conservation  of  mass  and  momentum  equations,  in  the  drift-diffusion  regime, 
coupled  with  the  Poisson  equation  for  the  self-consistent  calculation  of  the  electric  field.  An 
analytic  representation  of  the  force  density  on  the  neutral  gas  due  to  the  presence  of  a 
plasma  is  derived,  which  is  applicable  within  the  drift-diffusion  regime.  The  code  is 
validated  by  comparison  to  previous  calculations  and  then  applied  to  a  series  of  test 
problems  to  examine  the  behavior  of  the  plasma  force  density. 


D  =  Diffusivity  coefficient,  m2/s 
E  =  Electric  field,  V/m 

L  =  Time  rate  of  change  of  particle  density  due  to  loss  processes,  m~3  s'1 
P  =  Pressure,  ton¬ 
s'  =  Time  rate  of  change  of  particle  density  due  to  production  processes,  m'3  s'1 
T  =  Temperature,  K 

Particle  flux,  m"2  s'1 
Electric  potential,  V 
Permittivity  coefficient,  F/m 
Mobility  coefficient,  m2  V'1  s'1 
Frequency,  s'1 
Density,  kg/m3 


e  = 

Elementary  charge,  C 

r  = 

kB  - 

Boltzmann  constant,  J/K 

4>  = 

m  = 

Mass,  kg 

£  - 

n  = 

Number  density,  m'3 

P  = 

t  = 

Time,  s 

V  = 

v  = 

Velocity,  m/s 

p  = 

I.  Introduction 


Aerodynamic  flow  control  by  means  of  plasma  actuators  is  currently  an  active  area  of  research.  Such  devices 
have  experimentally  demonstrated  their  ability  to  reattach  separated  flow  at  high  angles  of  attack1’ 2’ 3  as  well  as 
to  induce  flow  movement  in  an  initially  stationary  air  mass4, 5.  The  use  of  plasma  actuators  in  such  a  role  may  offer 
several  advantages  over  traditional  flow  control  devices  (e.g.,  slats,  flaps,  slots).  Some  of  these  advantages  may  be 


-  Reduced  size  and  weight 

-  Increased  reliability6 

-  Increased  aerodynamic  agility 

-  Reduced  drag7 

-  High  bandwidth8 

-  No  moving  parts  (leading  to  improved  manufacturability)9 


*  Computational  Physicist,  Munitions  Directorate,  AFRL/MNAC,  101  W.  Eglin  Blvd,  Member  AIAA. 
f  Aeronautical  Engineer,  308  West  D  Ave.,  Member  AIAA 
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-  Inexpensive10 

-  No  tail  fins  (leading  to  increased  weapon  load  out  of  combat  aircraft)9 

Plasma  flow  control  devices  will,  no  doubt,  have  disadvantages  in  comparison  to  traditional  flow  control  devices 
as  well.  For  example,  plasma  actuators  may  not  perform  well  in  adverse  weather  conditions  and  will  require 
electrical  power  to  operate,  although  it  has  been  shown  that  this  latter  requirement  can  be  mitigated  through  a 
“smart”  actuator  methodology.11 

Whether  plasma  actuators  offer  any  demonstrable  advantages  over  traditional  flow  control  methodologies 
remains  to  be  seen.  However,  in  order  to  exploit  such  devices  to  the  fullest  extent  possible,  it  is  necessary  to  have  a 
basic  understanding  of  the  physical  phenomena  upon  which  they  are  based.  Although  the  experimental  component 
of  plasma  actuator  research  has  made  significant  discoveries,  the  modeling  and  simulation  of  these  devices  has  been 
conspicuously  lacking.  It  is  the  aim  of  the  present  paper  to  aid  in  this  respect  through  careful  computational 
modeling  of  the  relevant  physics  describing  the  heart  of  a  plasma  actuator,  the  dielectric  barrier  discharge  (DBD). 


II.  Physical  Model 

In  order  to  quantify  the  effect  of  the  plasma  on  the  neutral  flow  momentum,  a  suitable  momentum  transfer 
model  must  be  adopted.  The  model  proposed  in  this  work  is  given  as 


dA,v, 

dt 


+  V  p,  v, v,  =  -VP.  +  Y, P<v»  (v,  -  v„) 


(1) 


where  the  last  term  on  the  right  represents  the  time  rate  of  momentum  transfer  per  unit  volume  (i.e.,  force  density) 
from  electrons  and  ions  to  the  neutrals,  through  collisions.12  In  an  analogous  fashion  the  ion  momentum  equation 
can  be  modeled  as 


- + V  •  P,  v,  v,  =  -VP,  +  n,  e  E  -  p,  t/„  (v,  -  v. )  (2) 

where  the  second  term  on  the  right  describes  momentum  gained  by  the  ions  due  to  the  electric  field  and  the  last  term 
on  the  right  represents  momentum  lost  by  the  ions  due  to  collisions  with  the  neutrals.  A  similar  equation  can  be 
written  for  the  electrons.  The  simplest  model  describing  the  relevant  physics  of  a  DBD  is  based  on  the  drift- 
diffusion  model  of  a  plasma.  In  this  model,  the  formal  solution  of  Eq.  (2)  is  avoided  by  instead  solving  “pseudo- 
algebraic”  momentum  equations  for  each  species.  These  algebraic  equations  are  derived  by  neglecting  the  inertial 
terms  resulting  in  a  balance  between  the  electrical  force,  the  pressure  gradient  force,  and  the  collisional  momentum 
transfer  term.  Thus,  the  ion  velocity  may  be  approximated  as 


v,  =  v.  +• 


eE  VP, 


•",**  PiVi* 


(3) 


with  a  similar  result  for  the  electron  velocity.  Substituting  these  expressions  into  the  neutral  momentum  equation 
leads  to 


fov. 

dt 


+ V  •  p„  v«v»  =  -VP,  -  VPt  -  VP,  +  ( n ,  -nt)eE 


(4) 


fay. 

dt 


or 

+V-/JV  v  =  -VP 


(5) 
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with 


^=-'^Pt-^Pl+(n-n)eE  (6) 

Thus,  Eq.  (6)  defines  the  force  density  on  the  neutral  gas  due  to  the  presence  of  a  plasma.  The  presence  of  the 
plasma  adds  momentum  to  the  neutral  gas  in  three  ways:  via  an  electron  pressure  gradient,  via  an  ion  pressure 
gradient,  and  via  an  electric  field  term.  The  last  term  on  the  right  hand  side  of  Eq.  (6)  is  exactly  the  same  as  that 
given  by  Enloe,  et  al  (Ref.  4).  The  electron  and  ion  pressure  gradient  source  terms,  while  not  previously  reported  in 
the  plasma  actuator  literature,  can  be  shown  to  be  negligible  for  operating  conditions  that  are  typical  for  a  DBD  at 
atmospheric  pressure. 

The  rest  of  the  plasma  drift-diffusion  model  consists  of  the  equations  describing  conservation  of  mass  for  the 
electrons  and  ions 

^+V-(f^)  =  5eJ-4,  (7) 

and  Poisson’s  equation  which  describes  the  electric  potential 

V2<&  =  --(/2,-nJ  (8) 

e 

The  electric  field  and  the  potential  are  related  by  . 

£  =  -V<t>  (9) 
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distribution  is  small  compared  to  the  characteristic  time  for  the  discharge  to  develop.  This  is  considered  to  be  the 
case  for  DBDs  at  atmospheric  pressure  and  this  modeling  approach  has  been  widely  used.14,15,16'17,18’19 

In  summary,  the  plasma  drift-diffusion  model  consists  of  the  set  of  Eqs.  (7),  (8),  and  (9)  with  the  ion  and  electron 
fluxes  given  by  Eqs.  (10)  and  (13),  respectively.  In  principle,  the  solution  of  the  plasma  drift-diffusion  equations 
could  be  coupled  with  a  simultaneous  solution  of  the  neutral  gas  momentum  in  order  to  accurately  determine  the 
effect  of  the  plasma  on  the  neutral  flow.  In  practice  this  is  made  difficult  by  the  disparate  time  scales  corresponding 
to  the  plasma  actuator  and  the  neutral  flow.  For  example,  a  typical  plasma  actuator  has  a  length  of  approximately  1 
cm  and  induces  a  neutral  flow  velocity  of  approximately  3  m/s,20  resulting  in  a  characteristic  time  for  the  neutral  gas 
on  the  order  of  3  ms.  On  the  other  hand,  die  DBD  plasma  actuator  typically  operates  at  a  frequency  of  5-10  kHz, 
which  corresponds  to  a  characteristic  time  of  0.1 -0.2  ms,  or  a  factor  of  15-30  times  faster  than  the  neutral  flow.  This 
difference  in  time  scales  suggests  that  the  effect  of  the  plasma  on  the  neutral  gas  can  be  modeled  by  averaging  the 
body  force  term  in  Eq.  (5)  over  one  period  (7)  of  the  DBD: 

T 

|[-V^  (3c, t) -  VP,(x,  t)  +  (n,  (x,t)  -  nt  (5,  t))  e  E(x,  /)]  dt 
- - -  (i4) 

In  this  manner,  the  plasma  equations  are  decoupled  from  the  neutral  flow  equations,  and  the  effect  of  the  plasma  can 
be  accounted  for  by  an  average  body  force  term  in  the  neutral  gas  momentum  equation.  This  force  term  depends  on 
position  but  is  independent  of  time. 

At  this  point  it  must  be  mentioned  that  the  typical  plasma  actuator  is  inherently  at  least  a  two-dimensional 
device,  and  care  has  correspondingly  been  taken  in  writing  all  the  equations  presented  to  this  point  in  order  to  retain 
their  dimensionality.  However,  the  present  work  to  date  has  been  limited  to  a  single  dimension  to  examine  the 
physics  in  their  simplest  setting,  and  to  allow  for  a  systematic  numerical  building  block  approach.  Several 
researchers  have  pointed  out  the  difficulties  involved  in  capturing  the  correct  physics  with  a  one-dimensional 
approach,15,16,21  particularly  with  regard  to  the  calculation  of  the  electric  field.  Further  note  of  this  will  be  made 
where  appropriate. 


III.  Numerical  Scheme 

The  system  of  equations  (7,  8,  9)  is  highly  nonlinear  and  stiff,  and  therefore  numerically  difficult  to  solve.  The 
numerical  scheme  implemented  in  this  work  has  been  adapted  from  that  of  Boeuf,22  and  uses  an  implicit  difference 
representation  based  on  the  exponential  discretization  method  of  Scharfetter  and  Gummel.23  This  technique, 
originally  developed  for  simulations  of  electron  and  hole  transport  in  semiconductors,  is  robust,  stable,  and  able  to 
deal  with  situations  where  either  the  field  driven  flux  (first  term  in  Eq.  [13])  or  the  diffusive  flux  (second  term  in  Eq. 
[13])  is  dominant.  The  former  is  dominant  near  electrodes  while  the  latter  is  dominant  in  the  quasi-neutral  plasma 
region. 

The  temporal  integration  of  the  system  of  equations  is  performed  by  successively  solving  the  Poisson  equation 
(assuming  the  particles  do  not  move),  followed  by  solving  the  charge  continuity  equations  (assuming  the  electric 
field  does  not  change),  as  can  be  seen  in  Fig.  1.  Thus,  although  the  continuity  equations  are  solved  implicitly  in 
time,  the  overall  temporal  integration  scheme  is  explicit  in  nature. 

The  integration  of  the  continuity  equations  uses  a  simple  forward  difference  for  the  time  derivative  of  the 
number  density  and  a  central  implicit  difference  for  the  spatial  derivative  of  the  flux: 


ni  ni  |  A  14-1/2  1  M/2  _  gk 

A l  Ax 


(15) 


with  the  source  term  evaluated  explicitly,  as  shown. 
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Figure  1.  Overall  time  integration  methodology. 


The  Scharfetter/Gummel  exponential  spatial  discretization  is  used  to  represent  the  fluxes: 


(16) 


where 

07) 


where  s  is  equal  to  -1  for  negatively  charged  particles  and  +1  for  positively  charged  particles.  Note  that  the  number 
densities  are  evaluated  at  time  k  +  l,  while  all  other  quantities  are  evaluated  at  time  k.  Note  also  that  Z  is  unitless,  as 
it  should  be. 

The  main  advantage  of  this  scheme  is  that  it  provides  a  numerically  stable  representation  of  the  electron  and  ion 
flux  under  either  field-dominant  or  diffusion-dominant  conditions.  Representing  the  flux  using  a  standard  finite 
difference  approximation  leads  to  numerical  instabilities  when  the  electric  potential  between  adjacent  cell  nodes  is 
of  the  order  or  larger  than  the  characteristic  energy  (given  by  Einstein’s  relation  as  Z)///  =  kBT /e  ).22  On  the  other 
hand,  it  can  be  shown  that  use  of  the  present  method  reduces  the  discretized  PDE  of  Eq.  (15)  to  a  diffusion  equation 
(of  second-order  accuracy)  for  a  diffusion-dominated  flux,  and  to  a  convection  equation  (with  upwinding)  for  a 
field-dominated  flux. 

The  grid  is  defined  such  that  equally  spaced  nodes  are  located  throughout  the  interior  of  the  plasma  domain  (Fig. 
2),  with  cell  interfaces  located  halfway  between  cell  nodes.  Number  densities  and  electric  potential  are  defined  on 
the  cell  nodes,  while  the  fluxes  are  evaluated  at  the  cell  interfaces.  At  the  boundaries  of  the  plasma  domain,  the  grid 
is  modified  in  order  to  place  the  boundary  at  a  cell  interface  (Fig.  3).  This  methodology  has  useful  conservation 
properties24  and  makes  the  boundary  conditions  particularly  easy  to  implement. 


c* =tt  k+,A‘exP(z;ti/2)-cz)ii]  — ^ 


7* 

f^j+1/2 

exp(Z‘„2)-l 
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y-i/2 
-  1 


*+-1/2 


m - 1 - • 


Figure  2.  Discretization  of  the  computational  domain. 


Discretizing  the  fluxes  according  to  the  Scharfetter/Gummel  prescription  results  in  a  tridiagonal  system  of 
equations, 


Ak  *+!  ,  ryk  „t+l  ,  /~*k  *+l  _  n* 

nM+B,nt  +c,nM=Dl 


(18) 


which  is  efficiently  solved  in  0(3N)  operations  using  a  concisely  coded  procedure  based  on  LU  decomposition  with 
forward  and  back  substitution.25 


Boundary  conditions  must  be  applied  to  the  number  densities  and  fluxes  of  both  species  at  the  edges  of  the 
plasma  domain.  Accordingly,  the  normal  gradient  of  the  number  density  is  required  to  be  zero  at  a  boundary  for 
those  cases  in  which  the  electric  field  at  the  boundary  would  normally  drive  particles  out  of  the  plasma  domain. 
Therefore,  the  boundary  condition  for  outgoing  particle  flux  is  given  simply  by  Eq.  (10)  (reduced  to  a  field  driven 
flux  only,  since  Vn  =  0 ).  If  the  boundary  for  the  outgoing  particles  is  a  dielectric  surface,  the  flux  is  integrated  in 
time  in  order  to  calculate  the  surface  charge  density.  The  incoming  flux  for  positive  ions  is  set  to  zero  for  those 
cases  in  which  the  electric  field  at  the  boundary  points  into  the  plasma  domain.  The  incoming  flux  boundary 
condition  for  electrons  depends  upon  whether  the  boundary  is  an  electrode  or  a  dielectric  surface.  At  an  electrode 
boundary,  the  incoming  electron  flux  is  defined  to  be  a  fraction  of  the  outgoing  positive  ion  flux  at  the  same  location 
(but  opposite  in  sign),  i.e.. 


(19) 


where  y  is  defined  as  the  secondary  ionization  coefficient  (with  values  typically  between  0  and  0.5).  At  a  dielectric 
boundary  the  incoming  flux  is  described  by 


(20) 


where  is  the  desorption  frequency  for  the  electrons  leaving  the  dielectric  surface,  and  h  is  a  unit  vector 

pointing  into  the  plasma  domain.  This  description  is  similar  to  that  reported  by  Golubovskii.26 


3/2 


3 


Figure  3.  Discretization  of  the  computational  domain  near  a  boundary. 
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Poisson’s  equation  is  discretized  with  a  second-order  central  finite  difference  representation,  and  solved  using  a 
point  Gauss-Seidel  method  with  successive  over-relaxation.  In  practice  this  method  converges  in  a  few  thousand 
iterations  for  the  initial  transient,  thereafter  converging  in  only  a  few  iterations  per  time  step. 

For  cases  in  which  a  dielectric  barrier  is  present,  Laplace’s  equation  is  solved  in  the  interior  of  the  dielectric  and 
a  constraint  condition,  derived  from  Gauss's  law,  is  placed  on  the  normal  component  of  the  electric  field  at  the 
plasma/dielectric  interface: 


1  _  ^ ditltaric 

+  a 

dx 

plasma  ^ plasma  1 

£ 

dielectric  ** plasma 

(21) 


The  plus  sign  is  used  for  ID  configurations  in  which  the  dielectric  is  to  the  right  of  the  plasma,  and  the  minus  sign 
for  when  the  dielectric  is  to  the  left  of  the  plasma.  The  constraint  condition  is  included  as  part  of  an  iterative 
solution  for  the  potential  within  the  entire  spatial  domain.  That  is,  Poisson's  equation  and  Laplace's  equation  are 
alternately  solved  in  the  plasma  and  dielectric  domain  (respectively),  while  a  discretized  version  of  Eq.  (21)  is  used 
to  define  an  updated  value  of  the  potential  on  the  dielectric  interface  with  each  cycle.  This  process  is  repeated  until 
the  solution  through  the  entire  domain  converges. 

In  the  explicit  time  integration  approach,  in  which  successive  solutions  are  found  for  the  continuity  and 
Poisson’s  equation,  the  time  step  is  constrained  to  be  less  than  corresponding  dielectric  relaxation  time27-28  given  by 

t,  = - ^ -  (22) 

.  e( H'ti'+Hfr) 


IV.  Validation 

A.  Transient  Plasma  Sheath  for  a  Negative  Potential 

As  a  first  test  of  the  coupled  nlasma/potential  solver,  the  case  of  a  sheath  in  a  non-ionizing  Argon  plasma  at  a 
pressure  of  100  torr  is  examined.29  The  transport  parameters  used  are  as  given  as 

fit  =  0.3  m2/Vs  De  -0.3  m2/s 

/z,  =10"3  m2/Vs  £>.  =10"*  m2/s 

These  parameters  define  characteristic  temperatures,  based  on  the  Einstein  relation,  for  the  electrons  and  ions  of  1 
eV  and  0.1  eV,  respectively.  The  left  boundary  of  the  domain  is  an  electrode,  initially  at  a  potential  of  0  V.  The 
right  boundary  models  a  free  plasma  region  a  distance  of  200  characteristic  Debye  lengths  ( ZD0 ,  based  on  conditions 
in  the  unperturbed  bulk  plasma)  from  the  left  boundary.  The  electron  and  positive  ion  densities  are  spatially  uniform 
at  time  t-  0 ,  while  the  electric  potential  is  initially  zero  throughout  the  domain.  In  the  results  that  follow,  the 
various  quantities  have  been  nondimensionalized  by  their  characteristic  values:  number  densities  by  ne0 ,  electric 
potential  by  <b0  =  kDTt  !e ,  electric  field  by  £0=d>0/^D0,  the  force  density  by  F0=enc0E0,  and  time  by  the 
dielectric  relaxation  time  considering  only  the  ions  ( ta ).  At  time  t  =  0  the  left  electrode  is  instantaneously  brought 
to  a  nondimensional  potential  of  -50.  The  sudden  application  of  a  negative  potential  forces  the  highly  mobile 
electrons  away  from  the  electrode,  while  the  positive  ions  are  attracted  to  the  electrode.  As  the  charges  separate  an 
induced  electric  field  arises,  which  attempts  to  neutralize  the  total  field  (applied  +  induced)  within  the  bulk  of  the 
plasma. 
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(a)  Present  calculation  (b)  Previous  calculation  (Ref.  29) 

Figure  4.  Normalized  electron  and  ion  number  densities  as  a  function  of  distance  away  from  the  electrode  for 
test  case  1.  (b)  Results  from  a  previous  calculation  (Ref.  29)  under  the  assumption  that  the  electrons  are  in 
equilibrium  with  the  electric  potential,  (a)  Results  from  the  present  calculation,  which  makes  no  such 
assumption. 


A  comparison  of  number  densities  between  the  present  effort  and  previous  work  (Ref.  29)  is  highlighted  in  Fig. 
4,  which  shows  excellent  agreement  The  time  units  have  been  normalized  according  the  dielectric  relaxation  time 
corresponding  only  to  the  ions  (i.e.,  t/td,),  in  order  to  make  an  equal  comparison.  The  actual  calculation  was  made 


based  on  the  relaxation  time  given  by  Eq.  (22),  which  is  governed  by  the  electrons.  A  key  difference  between  Ref. 
29  and  the  present  effort  lies  in  the  treatment  of  the  electrons.  In  the  previous  work,  the  electron  continuity  equation 
was  not  solved;  instead  the  electrons  were  assumed  to  be  in  equilibrium  with  the  electrical  potential  throughout  the 
domain,  with  the  number  density  defined  according  to  the  Maxwell-Boltzmann  distribution: 


«.(*)  =  "«  expt 


(23) 


This  assumption  is  identical  to  neglecting  the  inertial  and  collision  terms  in  the  electron  momentum  equation,  and 
balancing  the  electrical  force  with  the  electron  pressure  gradient.  In  the  present  work,  this  assumption  is  not  made; 
both  the  ion  and  electron  continuity  equations  are  solved  simultaneously.  A  major  consequence  of  this  is  that  at 
t  ~  0  (just  after  the  normalized  potential  at  the  left  electrode  becomes  -50,  but  before  the  electrons  and  ions  can 
redistribute  themselves  to  cancel  out  the  total  electrical  field  in  the  bulk  plasma)  a  uniform  electric  field  exists 
throughout  the  domain  (see  Fig.  5a  for  t/t^  =1 ).  This  field  acts  to  drive  all  electrons  to  the  right  and  all  ions  to  the 

left.  Thus,  electrons  at  the  right  edge  of  the  computational  domain  are  lost  to  the  bulk  plasma,  while  positive  ions  in 
the  bulk  plasma  are  added  to  the  computational  domain.  The  highly  mobile  electrons  quickly  redistribute 
themselves  in  such  a  way  as  to  cancel  out  the  electrical  field  in  the  bulk  plasma.  The  time  scale  for  this  process  is  on 
the  order  of  several  tens  of  dielectric  relaxation  times,  as  given  by  Eq.  (22). 
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(a)  Present  calculation  (b)  Previous  calculation  (Ref.  29) 

Figure  5.  Electric  potential  throughout  the  domain  for  test  case  1.  (b)  Results  from  a  previous  calculation 
(Ref.  29)  under  the  assumption  that  the  electrons  are  in  equilibrium  with  the  electric  potential,  (a)  Results 
from  the  present  calculation,  which  makes  no  such  assumption. 


The  charged  particle  densities,  electric  potential,  and  electric  field  are  shown  in  Fig.  6a  at  time  tlt^  =  100 .  Near 
the  cathode  ( x-  0)  the  electric  field  is  negative,  and  there  is  a  region  of  net  positive  charge.  The  corresponding 
components  of  the  plasma  force  density  on  the  neutral  gas  (i.e.,  Eq.  [6])  are  also  shown  (Fig.  6b).  In  this  case  of  a 
sheath  in  a  negative  potential,  (i.e.,  one  in  which  the  sheath  is  located  in  a  potential  that  is  lower  than  that  of  the  bulk 
plasma)  it  is  observed  that  all  three  components  of  the  plasma  force  are  negative  and  that  the  force  due  to  the 
electron  density  gradient  can  be  significant  at  this  low  value  of  potential.  At  higher  values  of  applied  potential,  the 
force  diie  to  the  density  gradients  becomes  negligibly  small  in  comparison  to  the  electric  force. 


Figure  6.  Results  from  a  plasma  sheath  in  a  negative  potential  at  100  nondimensional  time  units.  AH 
quantities  have  been  appropriately  nondimensionalized. 
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A  space-time  plot  of  the  electric  potential,  electric  field,  total  charge  density  and  total  force  corresponding  to  the 
case  of  a  sheath  in  a  negative  potential  are  shown  in  Fig.  7.  It  is  observed  that  the  width  of  the  sheath  near  the 
cathode  ( x  =  0 )  grows  over  time,  and  that  the  direction  of  the  plasma  force  remains  in  the  negative  x-direction.  The 
magnitude  of  the  plasma  force  density,  however,  reduces  over  time.  This  test  case  will  be  revisited  later  in  this 
paper  in  order  to  examine  the  effects  of  a  positive  sheath  and  an  oscillating  sheath  on  the  plasma  force  density. 

B.  DC  Glow  Discharge 

The  code  is  next  validated  by  calculating  the  plasma  conditions  in  a  DC  glow  discharge  in  N2  at  a  pressure  of  5 
torn  The  cathode  is  grounded,  while  the  potential  applied  to  the  anode  is  determined  by  solving  the  circuit  equation 

K^=^~IR  (24) 

with  the  external  voltage  ( )  is  set  to  2  kV  and  the  resistance  (/?)  taken  to  be  300  kfl  Secondary  ionization  at  the 

cathode  is  allowed  with  an  emission  coefficient  (j)  of  0.1.  The  plasma  transport  and  ionization  parameters  are  as 
given  in  Ref.  21. 


Tinre(trt)  Time(t/T) 


Figure  7.  A  space-time  plot  of  the  electric  potential,  electric  field,  total  charge  density,  and  total  force  density 
corresponding  to  the  case  of  a  sheath  in  a  negative  potential.  The  cathode  is  located  at  xl  =  0 ,  and  the 
effective  anode  is  at  x! * lD0  =  200 .  All  quantities  have  been  appropriately  non-dimensionalized. 
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The  charged  particle  number  densities,  electric  potential,  and  electric  field  are  shown  in  Fig.  8a.  A  large  region 
of  positive  space  charge  is  present  near  the  cathode,  which  has  a  thickness  of  approximately  2  mm.  The  electric 
field  reaches  a  maximum  in  this  sheath,  with  a  value  on  the  order  of  4  kV/cm.  The  ion  density  also  reaches  a 
maximum  in  the  cathode  sheath,  due  to  the  high  ionization  rate  present  there  resulting  from  the  combined  effects  of 
the  high  electric  field  and  the  presence  of  electrons  from  cathode  emission.  Between  the  cathode  and  the  anode 
there  is  a  large  region  of  quasi-neutrality,  in  which  the  electric  field  is  relatively  low  (typical  value  -100  kV/cm).  A 
negative  space  charge  sheath  is  located  near  the  anode.  These  characteristics  are  typical  of  DC  glow  discharges  in 
general,  and  the  present  simulation  is  in  good  agreement  with  Ref.  30  in  particular,  which  was  calculated  under  the 
same  conditions. 

The  force  density  resulting  from  the  DC  discharge  is  shown  in  Fig  8b.  The  force  density  is  confined  to  the 
sheath  regions,  with  the  force  near  the  cathode  (where  a  large  electric  field  and  charge  density  difference  are 
present)  much  larger  than  the  force  near  the  anode.  In  both  cathode  and  anode  sheath  regions,  the  total  force  density 
is  due  almost  entirely  to  the  electric  field  and  is  directed  towards  the  respective  electrode. 

C.  RF  Glow  Discharge 

A  third  test  case  of  the  computational  code  is  that  of  an  RF  discharge,  with  the  working  gas  being  Helium  at  a 
pressure  of  1  torr.22  The  purpose  of  this  test  case  is  to  ensure  the  proper  behavior  of  the  code  in  ionizing  gases, 
correct  treatment  of  applied  voltages  at  high  frequency,  and  proper  inclusion  of  the  boundary  conditions  at 
electrodes,  including  secondary  ionization. 

The  right  electrode  is  grounded  while  the  potential  on  the  left  electrode  oscillates  sinusoidally  with  a  maximum 
amplitude  of  500V  and  a  frequency  of  10  MHz.  The  distance  between  electrodes  is  4  cm  (Fig.  9). 


(a)  Plasma  conditions  (b)  Force  density 

Figure  8.  Results  from  a  DC  glow  discharge  in  N2  at  a  pressure  of  5  torr. 


Figure  9.  Geometry  for  the  RF  discharge  used  in  test  case  2. 
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Initially,  the  electron  and  ion  densities  are  uniform  and  set  to  1015  m'3  everywhere.  The  electric  potential  is  also 
uniform  and  set  to  0  V.  Electron  and  ion  transport  coefficients  for  Helium  are  as  given  by  Boeuf22  and  Ward31.  The 
electron  (and  ion)  source  term  is  also  given  by  Ref.  31.  The  solution  procedure  is  to  iteratively  solve  the  electron 
and  ion  continuity  equations,  followed  by  the  Poisson  equation.  This  procedure  is  repeated  until  a  stationary  state  is 
achieved  in  which  the  periodic  form  of  the  solution  repeats  with  every  RF  cycle.  In  practice,  the  solution  takes 
several  hundred  periods  to  converge  as  such. 

Converged  values  of  the  electron  and  ion  number  densities,  as  well  as  the  electric  field  are  shown  in  Fig.  10  at 
four  specific  points  within  the  RF  cycle.  Note  that  at  10  MHz,  the  ion  distribution  remains  essentially  stationary, 
with  changes  in  the  electron  distribution  being  limited  to  the  sheath  regions.  The  agreement  between  these  results 
and  that  of  Ref.  22  (Fig.  1)  is  good. 

As  a  further  check  of  the  simulation,  the  electron  and  ion  conduction  current  and  the  displacement  current  on  the 
left  electrode  are  followed  in  time  and  compared  to  those  of  Ref.  22.  The  comparison  is  shown  in  Fig.  11;  again, 
agreement  is  good. 


Position  (cm)  Position  (cm) 


Figure  10.  Electron  and  ion  density,  along  with  the  corresponding  values  of  the  electric  field,  at  four  separate 
points  in  the  RF  cycle  of  a  1  torr  Helium  discharge  at  10  MHz. 
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Time  -vajrioit  ions  of  cunonl  dansities  and  applied  potential  at  the  left  electrode 


(a)  Present  calculation  (b)  Previous  calculation  (Ref.  22) 

Figure  11.  Time  variation  of  the  electron  (Je),  positive  ion  (Jjon)>  displacement  (JD)  and  total  (Jtoi)  current 
densities  on  the  left  electrode  of  a  1  torr  Helium  RF  discharge  at  10  MHz.  The  present  work  is  shown  in  (a); 
computational  results  from  Ref.  22  in  (b). 


V.  Results  and  Discussion 

Having  derived  an  expression  for  the  plasma  force  density  and  ensuring  the  code  is  working  properly,  we  next 
apply  the  simulation  capability  to  a  series  of  test  problems  in  order  to  examine  the  behavior  of  the  plasma  force. 

A.  Transient  Sheath  for  a  Positive  Potential 

The  sheath  in  a  non-ionizing  Argon  plasma  at  a  pressure  of  100  torr  is  reexamined  for  the  case  in  which  the 
sheath  is  located  within  a  positive  potential  (i.e.,  one  in  which  the  sheath  is  located  in  a  potential  that  is  higher  than 
that  of  the  bulk  plasma).  The  Argon  transport  parameters  remain  the  same  as  those  used  previously.  The  left 
boundary  of  the  domain  is  an  electrode,  initially  at  a  potential  of  0  V.  The  right  boundary  models  a  free  plasma 
region  a  distance  of  200  Debye  lengths  (based  on  the  unperturbed  bulk  plasma)  from  the  left  boundary.  The  electron 
and  positive  ion  densities  are  initially  spatially  uniform,  while  the  electric  potential  is  initially  zero  throughout  the 
domain.  At  time  t  =  0  the  left  electrode  is  instantaneously  brought  to  a  nondimensional  potential  of 
e&/kBTe  =+50.  The  sudden  application  of  a  positive  potential  draws  the  highly  mobile  electrons  towards  the 
anode,  while  the  positive  ions  are  repulsed  away  from  the  anode.  The  charge  density,  potential  and  electric  field  at 
tit =  100  are  shown  in  Fig.  12a.  Near  the  anode  the  electric  field  is  large  and  positive  and  there  is  a  region  of  net 

negative  charge.  This  field  acts  to  drive  electrons  to  the  left  and  ions  to  the  right.  Thus,  ions  at  the  right  edge  of  the 
computational  domain  are  lost  to  the  bulk  plasma,  while  electrons  in  the  bulk  plasma  are  added  to  the  computational 
domain.  The  ions  are  not  nearly  as  mobile  as  the  electrons,  thus  the  formation  of  the  sheath  occurs  much  more 
slowly  than  it  did  for  the  case  of  a  sheath  in  a  negative  potential. 

The  corresponding  force  densities  for  the  same  moment  in  time  are  shown  in  Fig.  12b.  Note  that  all  components 
of  the  plasma  force  on  the  neutral  gas  are  in  the  negative  x-direction,  just  as  in  the  case  of  a  sheath  in  a  negative 
potential.  That  is,  although  the  sign  of  the  potential  applied  at  the  electrode  has  changed,  the  direction  of  plasma 
force  on  the  neutral  gas  has  not. 

Space-time  plots  of  the  electric  potential,  electric  field,  total  charge  density  and  total  force  corresponding  to  the 
case  of  a  sheath  in  a  positive  potential  are  shown  in  Fig.  13.  It  is  observed  that  the  width  of  the  sheath  near  the 
anode  {x! XD0  =0 )  grows  over  time,  and  that  the  direction  of  the  plasma  force  remains  in  the  negative  x-direction. 
The  plasma  force  is  concentrated  near  the  electrode,  just  as  for  the  sheath  in  a  negative  potential.  The  magnitude  of 
the  plasma  force  is  reduced  over  time,  illustrating  the  transient  nature  of  the  plasma  sheath. 
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Figure  13.  A  space-time  plot  of  the  electric  potential,  electric  field,  total  charge  density  and  total  force 
density  corresponding  to  the  case  of  a  sheath  in  a  positive  potential.  The  anode  is  located  at  xl  =  0  and 
the  effective  cathode  is  at  xl  AD0  =  200 .  All  quantities  have  been  appropriately  non-dimensionalized. 
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B.  Transient  Sheath  for  an  Oscillating  Potential 

Having  examined  the  behavior  of  the  plasma  sheath  in  both  a  negative  and  a  positive  potential,  it  is  interesting  to 
examine  the  case  in  which  the  applied  potential  on  the  electrode  varies  sinusoidally.  In  this  case,  the 
nondimensional  potential  on  the  electrode  is  defined  as 

0(0,7)  =  50  sm{lnvt)  (25) 

with  v  =  5  kHz ,  a  typical  operating  frequency  for  a  DBD.  Space-time  plots  of  the  electric  potential,  electric  field, 
total  charge  density,  and  total  force  density  are  shown  in  Fig.  14.  In  these  plots  time  has  been  nondimensionalized 
by  the  oscillation  period  length  (0.2  ms).  It  is  observed  that  while  the  potential,  electric  field,  and  total  charge 
density  alternate  between  being  positive  and  negative  with  each  half-cycle,  the  plasma  force  always  remains  of  the 
same  sign.  As  in  the  previous  cases,  the  plasma  force  is  concentrated  near  the  electrode,  and  rapidly  diminishes  in 
magnitude  with  time. 

Caution  must  be  extended  if  one  attempts  to  directly  apply  these  results  to  the  case  of  a  DBD.  There  are 
numerous  differences  between  this  simple  test  case  and  an  actual  DBD:  in  the  present  case  the  model  is  ID,  a 
dielectric  barrier  is  not  present,  only  the  spatial  region  near  one  electrode  is  modeled,  and  the  voltage  is  quite  low. 
Additionally,  the  right  hand  side  of  the  computational  domain  in  this  model  is  open  (i.e.,  charged  particles  are  free  to 
enter  or  leave  the  domain  according  to  the  local  electric  field),  and  a  truly  periodic  solution  is  not  observed. 
Nevertheless,  some  insight  can  be  gained  by  examining  the  plasma  behavior  in  this  simple  model.  Notice  that  (after 
the  initial  transient)  the  plasma  force  "pushes"  the  neutral  gas  twice  each  period;  once  with  a  small  "push"  and  again 
with  a  larger  "push".  The  larger  push  is  associated  with  regions  of  net  positive  charge,  indicating  that  it  is  the  heavy 
ions  that  are  primarily  responsible  for  transferring  the  force  of  the  plasma  to  the  neutral  gas. 

C.  Microdischarges  without  a  Dielectric  Barrier 

Electrical  breakdown  in  gases  at  atmospheric  pressure  normally  occurs  through  a  large  number  of  individual 
breakdown  channels,  known  as  microdischarges.32  Many  researchers  have  studied  these  microdischarges,  including 
Dhali  (et  al,  Ref.  21),  who  examined  their  properties  in  the  absence  of  a  dielectric  barrier.  Their  two-dimensional 
simulation  is  repeated  here  in  one-dimension,  for  the  case  of  a  microdischarge  in  N2  at  atmospheric  pressure.  The 
physical  processes  of  ionization  (a),  recombination  (fi ),  diffusion  (D)  and  drift  (mobility  //)  are  included,  with 
corresponding  coefficients  given  as 

a  =  5.1  Pexp(-260P/ E)  cm-1  (with  pressure  P  given  in  torr  and  electric  field  E  given  in  V/cm) 

J3  =  lxl0-,A  m3/s 

=3.8158xl0'2  m2W  Dt  =1. 8000x1 0"'  mVs 

ju,  =  3.4211x10^  m2W  D,  =  2.3815xl(T5  m2/s 

Dhali  simulated  the  physical  process  of  photoionization  by  placing  a  homogenous  density  of  neutral  background 
ionization  throughout  the  domain.  This  was  not  performed  in  the  present  work. 

The  device  to  be  simulated  consists  of  two  plane  parallel  electrodes  a  distance  of  5  mm  apart,  with  the  left 
electrode  grounded  and  the  right  electrode  held  at  a  potential  of  26  kV.  As  an  initial  condition  the  electrons  and  ions 
are  Gaussian  distributed,  centered  on  the  cathode  with  a  1/e  point  of  0.27  mm.  The  electric  field  is  initially 
uniform  and  equal  to  -52  kV/cm. 

During  the  first  few  nanoseconds  an  anode-directed  streamer  is  observed  in  which  the  charged  particle  growth  is 
governed  by  Townsend's  first  ionization  coefficient  (Fig.  15a).  In  the  case  of  the  anode-directed  streamer,  the 
normal  drift  of  the  electrons  is  such  that  they  naturally  move  in  the  same  direction  as  the  streamer,  while  the  positive 
ions  move  in  the  opposite  direction.  This  results  in  a  net  negative  charge  at  the  streamer  head,  as  illustrated  in  Fig. 
16a,  taken  at  2  ns.  This  net  charge  density,  in  effect,  shifts  the  cathode  closer  to  the  anode,  increasing  the  electric 
field  near  the  streamer  head.  Thus,  the  streamer  propagates  into  an  ever-decreasing  electric  field,  which  serves  to 
both  increase  the  ionization  at  the  head  as  well  as  the  velocity,  as  evident  in  Fig.  15a. 
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Figure  14.  A  space-time  plot  of  the  electric  potential,  electric  field,  total  charge  density  and  total  force  density 
corresponding  to  the  case  of  a  sheath  in  an  oscillating  potential.  All  quantities  have  been  appropriately  non- 
dimensionalized. 


The  streamer  current  (Fig.  15b)  initially  decreases  as  the  electrons  and  ions  in  the  initial  Gaussian  distribution 
redistribute  themselves  in  order  to  cancel  out  the  applied  electric  field.  The  time  scale  for  this  process  is  of  the  order 
of  a  few  tens  of  dielectric  relaxation  times,  as  given  by  Eq.  (22).  The  exponential  rise  in  the  current  at  later  times  is 
indicative  of  Townsend  ionization  in  the  main  streamer  body. 

This  streamer  crosses  the  discharge  gap  in  approximately  5  ns,  resulting  in  an  average  speed  on  the  order  of  9.8  x 
105  m/s.  If  we  take  the  average  value  of  the  electric  field  at  the  streamer  head  as  it  crosses  the  gap  to  be  -78  kV/cm 
(a  value  determined  by  tracking  the  calculation  in  time),  then  the  average  electron  drift  velocity  at  the  streamer  head 
will  be  approximately  3.0  x  10?  m/s,  or  more  than  three  times  slower  than  the  streamer  velocity. 

Space-time  plots  of  the  electric  potential  (kV),  electric  field  (kV/cm),  total  charge  density  (cm'3),  and  force 
density  (N/cm3)  corresponding  to  a  microdischarge  in  N2  at  atmospheric  pressure  are  shown  in  Fig.  17.  The  non¬ 
zero  force  density  is  almost  entirely  localized  near  the  streamer  head,  and  in  contrast  to  previous  cases,  is  positive 
throughout  the  domain.  That  is,  the  force  is  continuously  directed  towards  the  anode,  even  when  the  streamer  is  far 
from  the  anode.  This  is  due  to  the  large  negative  electric  field  that  exists  at  the  streamer  head,  combined  with  the 
net  negative  space  charge  that  is  also  present  there.  The  positivity  of  the  force  indicates  that,  on  a  time  scale 
characteristic  of  a  microdischarge,  electron-neutral  collisions  are  the  dominant  mechanism  for  transferring 
momentum  from  the  plasma  to  the  neutral  gas.  For  example,  in  the  time  it  takes  for  the  microdischarge  to  traverse 
the  entire  gas  gap  (5  mm),  the  heavy  ions,  with  their  much  lower  mobility,  will  have  moved  a  distance  of  the  order 
of  0.02  mm,  at  most  (assuming  an  average  electric  field  of  -78  kV/cm,  as  determined  previously).  The  actual 
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Figure  15.  Profiles  of  electron  density  (a)  and  current  (b)  for  a  microdischarge  in  N2  at  atmospheric  pressure. 
The  applied  voltage  on  the  anode  (located  at  x  =  5  mm)  is  fixed  at  26  kV  and  the  initial  charged  particle 
densities  are  Gaussian  distributed,  centered  on  the  cathode  with  a  1/e  point  of  0.27  mm. 

distance  traveled  by  the  ions  during  this  time  will  be  much  less  than  this  because,  unlike  the  electrons  at  the  streamer 
head  which  naturally  move  in  a  direction  of  a  greater  electric  field,  the  ions  will  move  toward  the  "back  side”  of  the 
streamer,  which  is  a  region  of  much  reduced  (in  magnitude)  electric  field  (see  Fig.  16). 


(a)  t  =  2  ns 


(b)  t  =  4  ns 


Figure  16.  Profiles  of  particle  density,  potential  and  electric  field  for  a  microdischarge  in  N2  at  atmospheric 
pressure  at  (a)  2  ns  and  (b)  4  ns.  These  times  correspond  to  an  anode-directed  streamer  and  a  cathode- 
directed  streamer,  respectively.  The  applied  voltage  on  the  anode  (located  at  x  =  5  mm)  is  fixed  at  68  kV  and 
the  initial  charged  particle  densities  are  Gaussian  distributed,  centered  on  the  cathode  with  a  Me  point  of 
0.27  mm. 
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Figure  17.  Space-time  plots  of  the  electric  potential  (kV),  electric  field  (kV/cm),  total  charge  density  (cm"3) 
and,  total  force  density  (N/m3)  corresponding  to  a  microdischarge  in  N2  at  atmospheric  pressure.  The 
cathode  is  at  x  =  0  mm,  the  anode  at  x  =  5  mm,  and  the  applied  potential  is  68  kV. 


VI.  Conclusion 

Several  simulations  have  been  made  of  various  plasma  devices.  These  simulations  are  based  on  successive 
numerical  solutions  of  the  electron  and  ion  charge  conservation  equations  (in  the  drift-diffusion  approximation) 
coupled  with  Poisson’s  equation.  The  code  was  validated  with  three  test  problems,  illustrating  its  ability  to  correctly 
solve  the  governing  equations,  and  then  applied  to  several  simple  problems  to  examine  the  behavior  of  the  resulting 
force  density.  Future  work  will  consist  of  extending  the  model  to  higher  dimensions,  examining  the  effect  of 
negative  ions  (present  in  air  plasmas),  and  optimization  studies  involving  electrode  geometry,  dielectric  permittivity 
and  applied  voltage  waveform. 
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